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ABSTRACT 


The objective of this research was to develop an optimization technique that can be 
used interactively by design engineers to approach an optimal! design with minimal 
computational effort. The technique can be applied to both continuous and discrete values 
of design variables. A large number of design variables can be also considered. 

In order to meet the objective, an optimization procedure was developed by coupling 
the finite element analysis (FEA) to the orthogonal array experimentation technique, because 
FEA is a common analysis tool for design engineers. From the results of the FEA and an 
orthogonal array, an average Jacobian matrix was constructed that showed the average 
overall sensitivity of the design variables. These sensitivities were then used to optimize the 
design parameters. The process could then be repeated at the discretion of the engineer until 
a satisfactory design is obtained. In general, the designer can predict and control the 
number of FEA calculations before an optimization process so that one can plan a budget 
and time for an optimal design. 

Some examples of structural optimization with truss and frame structures with 
continuous and discrete values of design variables were studied using the technique 
developed in this paper. Their optimal solutions were found with small numbers of 


iterations. 
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I. INTRODUCTION 


Optimization is a process that can be seen in almost every aspect of life. Humans are 
always striving to go faster, climb higher, run farther, and engineering is no different. 
Engineering requires optimization so that we can design faster planes and cars, that use less 
fuel, are lighter, stronger, and more comfortable. Engineers have been optimizing designs 
since the beginning, but recent advances in computing have made numerical optimization 
techniques a more effective way than the original trial-and-error and experience based 
optimization. The present methods based on numerical search techniques are effective and 
easy to understand and implement on design problems that don't have many design 
variables (Vanderplaats, 1984). However, it is difficult to estimate the computational time 
and budget to converge to an optimal design using most of the numerical optimization 
techniques available today unless the design engineer has experience with the same kind of 
design problem in advance. The computational costs increase highly nonlinearly as the 
number of design variables increases. Therefore, application of a numerical optimization 
technique to a practical engineering design has been hampered. 

The objective of this research was to develop an efficient method for achieving the 
optimal design with minimal computational effort. Also, the developed technique allows the 
engineer to interact with the procedure and to estimate how much time is to be spent on the 
optimization depending on the number of design variables (it does not matter how many 
design variables there are), and how close to the actual optimum the design must be, 
balancing schedule and cost. The approach was to develop an optimization procedure by 
combining the finite element method and orthogonal arrays. The emphasis was in structural 
mechanics, but could be applied to other disciplines. Because finite element analyses are 


time consuming in general, the approach was to minimize the number of iterations 





associated with the FEA. This led to the incorporation of orthogonal arrays in the 
methodology in order to reduce the number of FEA computations. 

Orthogonal arrays were developed initially for discrete values of design variables. 
Therefore, they work well for the discrete case. Further, the optimization technique was 
also applied to continuous design variables with design constraints. 

The goal is not necessarily to get the actual optimum value, assuming that one exists, 
but to develop a more practical approach that will achieve a result close to the optimum with 
reasonable effort. In addition to developing a methodology that would lead to an optimal 
design, it was necessary to identify the overall sensitivity of the various parameters. This 
information provides an engineer with further physical understanding of the design 
problem. 

Although the technique may not achieve the most optimal value, it will give the 
designer the value of parameters close to the optimum with reasonable effort and drastically 
reduced computational time. Just like if you want a higher quality car you might buy a 
luxury car, or a faster car you might buy a sports car, you must pay the price. Similarly, if 
you want the most optimum design you must pay for it with either higher design and 
manufacturing costs or more computer time. On the other hand, if you are willing to accept 
a good solid car you might buy a basic family car. This is similar to spending a reasonable 
amount of time optimizing and getting a close to optimal design without the high computer 
computational times. The latter approach is therefore the goal which should lead to a more 
practical optimization approach. | 

This paper will first discuss the background of the Taguchi Method and the 
orthogonal arrays that he has incorporated into his design of experiments methodology. 
This will be followed by a brief technical description of the orthogonal arrays as described 
by Peace (1993). Next will be the method and procedure developed in this research. The 


final section will include some examples of the technique. 





II. TAGUCHI METHOD AND ORTHOGONAL ARRAYS 


The optimization approach developed in this research centered around the concepts 
developed by Taguchi for quality engineering and then coupling his concepts with the finite 
element method to achieve an optimal design methodology. The following brief description 
of orthogonal arrays can be found in Peace (Peace, 1993). 

The orthogonal array was the foundation of the methodology that Taguchi developed 
for designing an experiment. Although other methods of conducting an experiment are 
possible, like a full factorial approach, the orthogonal array has traditionally been associated 
with the Taguchi experimentation technique. The orthogonal array is repeatable, cost 
effective, and efficient in that only a small amount of data must be collected and then that 
data can easily be translated into meaningful and verifiable conclusions (Peace, 1993). 

The orthogonal array was originally developed by Fisher of England. It was 
developed during his efforts to control error in an experiment (Ealey, 1988). Later, Taguchi 
adapted the orthogonal array in order to measure the effect of a factor being studied on the 
average result (Peace, 1993). 

The details of the orthogonal array will follow, but the concept of orthogonality 
refers to the statistically independent or balanced parameters that make up the columns of the 
orthogonal array. One major objective of this research was the application of the technique 
in both discrete and continuous parameter problems. However, the orthogonal array 
approach lends itself directly only to problems with discrete parameters. For problems 
associated with discrete parameters, the orthogonal array is used directly. But, for problems 
with continuous variables, this research will show that taking the minimum and maximum 
values in the range of each of the parameters, and then putting those values into the 
orthogonal array, is an effective way of dealing with continuous variables. A discussion of 


how this was dealt with will follow this section. 





In a typical orthogonal array as shown in Table 1 (Taguchi, 1987), each level (here, 
1 or 2) has an equal number of occurrences in each column. Also notice the relationship 
between one column and another. For each level within one column, each level within any 
other column will occur an equal number of times as well. This introduces the statistical 
independence or balance into the orthogonal array. Since each column is orthogonal to the 
others, if the results associated with one level of a specific factor are much different at 
another level, it is because changing that factor from one level to another has a strong impact 
on the design parameter being measured. Because the levels of the other factors are 
occurring an equal number of times, any effect by these other factors should be canceled out 


(Peace, 1993). 


Table 1. Lg(2”) Orthogonal Array 
Ls(27) 





Parameter No. | 1 






No. of Runs 
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The design of an orthogonal array experiment, although balanced, does not require 
that all combinations of all factors be tested. The full factorial experiment studying 7 factors 
would require 128 (= 2’) experimental runs, compared with the 8 runs in the orthogonal 
array. The advantage of the orthogonal array can be seen in its simplicity, time savings, and 


therefore cost savings. 





It should be noted, that the columns in the orthogonal array can represent not only 
design factors, but also interactions. Therefore, by using sound engineering judgment, the 
orthogonal array can include critical interactions. The orthogonal array, does not study all 
interaction combinations like a full factorial experiment, but the balanced nature of the 
orthogonal array can compensate for this. 

The terminology of the orthogonal array La(B°) is as follows. The subscript of L, 
which here is A, represents the number of experimental runs which can be conducted in the 
experiment. B denotes the number of levels within each column, and C denotes the number 
of columns (factors and/or interactions that can be studied) in the orthogonal array. 
Therefore, the above orthogonal array Lg(27) contains eight experimental runs. Each 
column contains '2' levels, meaning each factor will posses two levels. The '7' means that 


up to seven factors and interactions can be incorporated into the experiment. 

















Ill. METHOD & PROCEDURE 


As discussed earlier, the procedure that has been developed focuses on the 
orthogonal array. The following example is written for a structural optimization problem 
with stress, displacement, natural frequencies, etc. as the design constraints, but in general 
could be applied to other disciplines. The design variables varied from a minimum value to 
a maximum value continuously or discretely. The objective was to reduce the weight of the 
structure supporting a load while not exceeding the design constraints in any of the members 
of the structure. In this example, the cross-sectional areas are considered as the design 
variables, the weight of the structure was the objective function, and the stresses, 
displacements, and natural frequencies were the constraints. 

The general procedure of the optimization technique that was developed begins by 
constructing the proper orthogonal array for the given design variables using the entire range 
of values for each variable. If the design variables are continuous, then the values of 1 and 
2 are used in the orthogonal array corresponding to the minimum and maximum values 
respectively. If the design variables are discrete, then the lowest discrete value is used for 
level 1 in the orthogonal array, the next lowest value is used for level 2 in the orthogonal 
array, etc. 

The next step is to use the design values in each row of the orthogonal array and run 
a FEM program (Kwon, 1997) for each row using the corresponding design values in the 
columns of the orthogonal array. The inputs to the FEM code are the cross-sectional areas, 


and the output from the FEM code are the values of normalized stress, displacements, 


natural frequencies, etc., P,. Here P. is the constraint function and is defined as 


p =Sticy (1) 








where f; is a computed value from the finite element analysis such as stress, displacement, 
or natural frequency. And, S, is the allowable value of f, such as failure strength, 
displacement, or natural frequency, and k, is a factor of safety (k, 21). 

After all cases have been run in the orthogonal array, the results are plotted P, vs. x j 


which is the sensitivity of a normalized constraint value with respect to a design variable, 


namely, the average overall sensitivities (slopes) ra are calculated. With the desired value 
j 


for all P, = 1, a target value of P. was selected using a one-dimensional (1-D) search 
technique in the direction determined by the average overall sensitivities until a constraint 
boundary was encountered. 

The next step is to take the best value from all of the orthogonal array runs that meet 
the design criteria and constraints, although it will be shown that the final solution is 
relatively independent of the starting point. It is recommended that the starting point be in 
the feasible domain for all parameters, because the technique is attempting to approach the 
optimum design from within the feasible domain. If no cases in the orthogonal array results 
are feasible, the most conservative design may be used as the initial design. These cross- 
sectional areas will then form the vector, {x/}. These values are then used to run the FEA 
code, which produces a vector of P-values called {PJ}. The superscript denotes the iteration 
number. 

From the orthogonal array, the average overall sensitivity of each parameter over its 
entire domain is calculated. This is slightly different from the common optimization 
approach of calculating the Jacobian (which contains the local sensitivities calculated at one 
point). In the traditional approach, the Jacobian matrix is determined with a fixed value of 


each design variable. The Jacobian matrix (Lang, 1968) is 





eo ee 
dx, dx, dx, 
gh > ., Ma 
[J]=| dx, ax, a, (2) 
dP, dP, _ bP, 


where n is the number of design variables and m is the number of constraints. In the 
present study, the averaged Jacobian matrix over the domain of each variable can be 
assembled from the average overall sensitivities instead of being computed at a specific set 


of design values. Thus the averaged Jacobian matrix is 

















AP, AP, AF, 
Ax, Ax, Ax,, 
AP, AP, AP, 
[J]=| Ax, Ax, Ax, (3) 
AP. AP. AP. 
Ax, Ax, Ax, 


In general , the off-diagonal terms in the first n rows are small compared with the 


diagonal terms and are therefore neglected. Thus [J | can be rewritten as 





























APL og 0 
Ax, 
0 sale 0 
Ax, 
- : ‘ AP. 
[J | =| 0 Ax. 
AP, 41 AP 41 AP. 41 (4) 
Ax, Ax, Ax, 
AP,, AP, AP, 
Ax, Ax, Ax, 


Out of the terms in the averaged Jacobian matrix, the absolute maximum in each column is 
used because that value is corresponds to the most constraining value when the design 
variables are updated. Using these sensitivities, the following equation is solved to compute 


the change of the design variables 
a7, (22) . 65) 


for i=1...n and j=1...m. Here Ax; is simply the needed change in each design variable to 


achieve an optimal value which is to be solved for using 


—] 
AP. 
Ax. = AP.| —+ 6 
| (2) 6 
where 
AP. = {P, — P,} (7) 
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Here, AP. is the distance to be traveled in the direction determined by the orthogonal array 


from the current value of P. to the target value P.. The target value P. is the value to be 


determined from the 1-D search. This value can be found by extrapolation until a constraint 
is violated and then interpolated until the constraint boundary is found. Many 1-D search 
techniques exist that can efficiently accomplish this task such as the bisection method or the 
golden section method (Vanderplaats, 1984). 

Ax; is then added to the previous design value x;! to yield the new design values X;*. 
The new design values are then used in running the FEA program again to verify the result. 

At this point, if any or all of the values in P;? are close to 1, then this is assumed to 
be very close to the optimum design, and the analysis can be stopped, with the design 
values becoming x;*. If further reduction is required, then the range of the parameter values 
is refined by using +5% of the values in x;. The procedure is then repeated using the new 
refined maximum and minimum values as the limits in the orthogonal array until a 
satisfactory design is obtained. Usually, more than one set of orthogonal arrays is not 
necessary unless a very refined optimal design is required. This further refinement is the 


topic of follow-on studies. 


For a discrete value of a design variable, the change of design variable Ax; is 


compared to the allowable change of the design variable. If the magnitude of the calculated 


Ax; is greater than the magnitude of the allowable change, then Ax; is set equal to the 
allowable value. Otherwise, the design variable does not change. This process continues 
until there is no longer any change in the design. 


Figure 1 outlines the procedure described above. 


1] 





Optimization Procedure using Orthogonal Arrays 


. Construct the Proper Orthogonal Array. 
. Select an Initial (Feasible) Design. 


. Calculate the Average Overall Sensitivities. 


. Update the Design Variables using a 1-D Search Technique. 


. Ifthe design is satisfactory, then stop the procedure. Otherwise, repeat 





the procedure with a refined domain for further optimization. 


Figure 1. Procedure Summary 
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IV. RESULTS & DISCUSSION 


In this section, three examples will be shown to illustrate the technique described in 
the previous section. The first example is a statically independent 3-bar truss with 
continuous design variables and stress/displacement constraints. The second example is a 
statically independent 6-bar truss in a star shape with discrete design variables and only 
stress constraints. The third example is a 12-bar frame structure with continuous design 
variables. It has stress, maximum displacement, and minimum natural frequency 
constraints. The frame problem in the third example is first solved using a conventional 
numerical optimization technique, Sequential Quadratic Programming (SQP). Next, the 
frame problem is solved with orthogonal arrays using the techniques developed in this 


research. Finally, the results of the two optimization techniques are compared. 


A. STATICALLY INDETERMINATE THREE-BAR TRUSS EXAMPLE 





Figure 2. 3-Bar Truss 


This first example is a statically independent three-bar truss structure subjected to a 
load of 50,000 Ibs to the right and 100,000 lbs down. The objective is to minimize the 


weight of the structure such that each member does not exceed the failure strength of 118 ksi 
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and the displacement at the node where the load is applied should not exceed 0.12 inches 
horizontally and -0.06 inches vertically. A factor of safety of 1.5 was used for the stress 
and displacement values. Initially, the minimum cross-sectional area is 0.25 in2, and the 
maximum cross-sectional area is 2.00 in? for all members. The design variables (cross- 
sectional area) in this problem were continuous. These areas are then put into the following 


L4(23) orthogonal array (Phadke, 1989) (See Table 2). 


Table 2. L4(23) Orthogonal Arra 


Area ] Area 2 Area 3 
Run # 


1 0.25 0.25 0.25 
0.25 2.00 2.00 





























0.25 2.00 
2.00 0.25 


2.00 
2.00 





2 
3 
4 





Table 3 shows the results from running the FEM code for each of the 4 rows in the 


orthogonal array, Table 2. 


Table 3. FEA Results 
Pl 
1.03 1.69 4.33 
0.27 0.32 0.49 





2 
3 0.28 0.35 0.62 
4 -0.15 0.52 3.22 


The first three columns (P1 - P3) are the normalized stresses and the next two 
columns (P4 - P5) are the normalized displacements of the node where the load is applied. 

Figure 3 is a sensitivity plot for the normalized stress in element 1 with respect to the 
three different cross-sectional areas (Al, A2, A3). Similarly, Figure 4 is a sensitivity plot 


for the normalized stress in element 2 with respect to the three different cross-sectional 
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areas. And finally, Figure 5 is a sensitivity plot for the normalized stress in element 3 with 
respect to the three different cross-sectional areas. In all these three plots, the dotted lines 
represent the stress relative to element 1, the dashed lines represent the stress relative to 
element 2, and the dash-dot lines represent the stress relative to element 3. From these 
plots, one can see that the slope of the normalized stress in each member vs. the area of that 
same member is the steepest (or close to the steepest) and is therefore has the highest 
average overall sensitivity. 

We then take only the sensitivity plot for element 1 with respect to the cross- 
sectional area of element 1 and plot this sensitivity in Figure 6. Similarly, Figure 7 is the 
sensitivity of element 2 with respect to the cross-sectional area in element 2. Finally, Figure 
8 is for element 3. Notice that the slope in Figure 8 is the steeper than Figures 6 and 7, 
indicating that the stress level in element 3 is much more sensitive to changes in cross- 
sectional area 3 than in elements 1 and 2. 

Figures 9 - 11 show the sensitivities of the normalized horizontal displacement, P4, 
to the three cross-sectional areas. Here we see that all of the slopes are about the same 
indicating that the horizontal displacement is affected equally by changes in all three 
elements. 

Figures 12 - 14 show the sensitivities of the normalized horizontal displacement, P4, 
to the three cross-sectional areas. Figure 14 shows the steepest slope, indicating that 
changes in the cross-sectional area of element 3 will have the greatest impact on the vertical 


displacement. 
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Pi vs. Area 


P = Stress / Strength 
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Figure 3. Sensitivity Plot, Pl vs. Area 


P3 vs. Area 


P = Stress / Strength 
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Figure 5. Sensitivity Plot, P3 vs. Area 
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P2 vs, Area 
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Figure 4. Sensitivity Plot, P2 vs. Area 
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Figure 6. Sensitivity Plot, Pl vs. Al Figure 7. Sensitivity Plot, P2 vs. A2 


P3 vs. A3 
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Figure 8. Sensitivity Plot, P3 vs. A3 
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P4 vs. Al ene 
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Figure 9. Sensitivity Plot, P4 vs. Al Figure 10. Sensitivity Plot, P4 vs. A2 


P4 vs. A3 
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Figure 11. Sensitivity Plot, P4 vs. A3 
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P5 vs. A2 


P5 vs. At 
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Figure 12. Sensitivity Plot, P5 vs. Al Figure 13. Sensitivity Plot, P5 vs. A2 


P5 vs. A3 





P5 = V/ Vmax 
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Figure 14. Sensitivity Plot, PS vs. A3 
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In this analysis we will use run 3 as the starting point because it has the lowest 


weight of the runs where no constraints are violated. The starting point {A!} is 


{a'l’ ={2.00 0.25 2.00} (8) 
with the resulting 
{P'\" ={0,.28 0.35 0.62 0.41 0.91} (9) 


It meets the criteria that none of the constraints are violated, therefore the design variables 
are all in the feasible domain. This run results in a weight of 78.71 Ibs. This becomes the 
{A!}vector with the corresponding result, {P!}. 

In order to optimize the design variables, the following calculations were performed 
for each member. For example, for member 1 equation (5) is applied for both stress and 
displacement constraints with a target value of 0.9. From this, three values (1 from stress . 


constraint and 2 from the displacement constraints) of AA were calculated, then the largest 
absolute value was chosen to update the design variables. The resulting design variables 


were 


{a?\" ={0.75 0.25 1.90} (10) 


In order to verify that all of the design constraints were met, the FEM code was mun using 


these design variables, with the following results 


{P?}' ={0.65 0.65 0.63 0.97 0.94} (11) 
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The corresponding weight for this design is 47.55 Ibs, a considerable weight reduction after 
only five FEA computations. The results show that the new design variables satisfy the 
constraints and are therefore an acceptable solution. However, if the engineer wants to 


further optimize the design variables, then the procedure can be repeated with a refined 


domain. 


The refining of the domain is done individually for each cross-sectional area. 


{ Ayan} is calculated from the equation 
{Ali} = {4/"}- B{(1-|P)-{477}} (12) 
similarly, {Ayax} is calculated from the equation 

{Ale} ={4""} + B{(1-|P*")-{4’"}} (13) 


where B is a scaling factor to adjust the refined domain. Once again, the upper and lower 
values are checked to insure that the values remain in the feasible domain. 

In general, the procedure using the orthogonal array only needs to be carried out 
once in most practical applications. This results in a reasonably optimum design without 
exhaustive computations. The final results of this procedure are relatively independent of 


the starting point in the analysis. 
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B. STATICALLY INDETERMINATE STAR-SHAPED TRUSS EXAMPLE 


Figure 15. 6-Bar Star-Shaped Truss 


Following the procedure outlined in the previous section, this example shows a 
statically indeterminate six-bar truss structure with discrete design variables. The geometry 
was chosen to insure that the problem was fully statically indeterminate. The objective here 
was to minimize the weight of the structure, while not exceeding the failure stress in any of 
the six members. The length of all 6 members was 10 in. The discrete cross-sectional areas 
for all six members were 0.25 in? , 0.50 in?, 1.00 in?, 2.00 in2, and 2.50 in2. These 
represent the five discrete design variables to be used in the orthogonal array. The L5(59) 
orthogonal array (Taguchi, 1987) was used for this problem. Table 4 shows the FEA 


results for each run. 
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Table 4. FEA Results 
Pl 


1.11 
1.41 
0.96 
1.12 
0.96 
0.79 
1.04 
0.85 
0.67 
0.95 
0.72 
0.61 
0.85 
0.91 
0.72 





From the results shown in Table 4, the overall sensitivities of the constraints to the 





design variables — between each level, are computed. For each level apply Equation 5 
Jj 


and then compute AA; . If AA; is less than the amount available for the reduction of each 
design variable, then the value of the design variable can be further reduced. 

In this example, the 7th run is used as the initial design and is updated based upon 
the previously described procedure. Notice that the weight of the 7th run is 67.5 Ibs., 
which ends up being very close to the final value. One could stop here and get a very good 


feasible design. The final result was 
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fa?\" ={0.50 0.25 1.00 2.00 2.00 0.25} (14) 


{p*\" ={0.97 -0.58 0.62 0.92 0.76 -0.16} (15) 


Because each of the values in the results vector does not violate any of the design 
constraints, the values in Equation 14 are acceptable design variables. The weight for this 
structure is 60 Ibs., a reduction of 40% over the most conservative design for this structure 


of 150 Ibs in 26 FEM calculations. 


C. OPTIMIZATION OF A 12-BAR FRAME STRUCTURE 


As discussed earlier, the 12-bar frame structure problem is first solved using 
common numerical optimization techniques. Next, the problem is solved using the 


procedures developed in this research. Figure 16 shows the 12-bar frame structure. 
1. Solution Using Common Numerical Optimization Techniques 


The purpose of this example was to minimize the volume of a frame structure 
subject to a horizontal load of 1000 Ibs at the top of the left edge, and 1000 Ibs. at the 
center of the left edge (See Figure 16). The constraints on the problem included the 
maximum stress in any member could not exceed the yield strength of the material, the 
maximum displacement of any node was not to exceed .125", and the lowest natural 


frequency of the system must be larger than 7 Hz. 
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Figure 16. 12-Element Frame Structure 


The cross-sectional area of each member was rectangular with a 2:1 height-to-width 
ratio. Therefore, the heights of each of the members made up the 12 design variables. The 
upper limit for each member was h=4.0", and the lower limit for each member was h=0.1". 

The main program executed the MATLAB constr optimization routine that called a 
‘function M-file' (fun.m) which calculated the values of the objective and constraints 
(Grace, 1992). The ‘function M-file’ also called 2 FEM programs that were used to 
calculate the constraints. The first FEM program (disp_frame_fem.m), calculated the stress 
in each member and determined the maximum horizontal displacement of any node. The 
second FEM program (eig_frame_fem.m), calculated the natural frequencies of the system. 

The Brobleih statement for the constrained optimization problem is given by the 


following objective function, F(h), and constraints, g(h): 


. 12 
minimize: F(h) = > bh, L (16) 


i=] 


2 























subjectto: gi) =|2™"|-1<0 fori = 1....,12 (17) 
yield 
@ 
13)=1-—*<0 18 
g(13) 0 (18) 
2(14)=—*+-1<0 (19) 
0.1<h, <4.0 for i = 1,...,12 (20) 
where: Umax = 0.125" 
(2 = 7 Hz 


Syield = 78000 psi 
bj = 0.5 h; 
h; are the 12 design variables 


The constraints are calculated from FEM calculations 


The constr function in MATLAB used for this constrained optimization problem, 
uses a Sequential Quadratic Programming (SQP) method. In this method, a Quadratic 
Programming subproblem is solved at each iteration, and an estimate of the Hessian matrix — 
of the Lagrangian is updated at each iteration using the BFGS formula. (Grace, 1992) The 
termination criterion used for this program are as follows: 

Termination criteria for h = 1x10! 

Termination criteria for f = 1x10"! 

Termination criteria for g = 1x10-2 


Maximum function calls = 1000 


The program was run for 6 different initial values of hp. Table 5 shows the value 


of the objective function at the optimum in each case. Also shown are the values of the 14 
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constraints and the 12 design variables corresponding to the optimum. The last row in the 
table shows the number of function evaluations for each trial. The function evaluation 
denotes a set of finite element analyses of the frame structure for both static and eigenvalue 
analyses. 

The results show that the design was improved tremendously over the initial 
design. The objective function for the 4.0 initial design (worst case) was, F(h®) = 720.00 
in3. Compared with the corresponding optimum design of F(h*) = 230.69 in?. 

From the results in Table 5, it can be seen that depending on the starting point, the 
optimum design goes towards 1 of 2 directions. Initial values of 4.0, 3.9, 1.0, and 0.1 
converge towards 1 final design. While, initial values of 3.0 and 2.0 converge towards 
another final design. The different designs can be characterized by the h(1) and h(12) 
values. Both designs have essentially the same value of the objective function, however 
the designs are considerably different. It appears that either the left side or the right side of 
the frame, but not both, carries the bulk of the load depending on the initial design. 

An additional constraint of the buckling load in each compressive member would be 
useful and might alleviate the problem of the 2 final designs. Also, constraints that restrict 
the relative displacement of the nodes might prove useful. Further improvement in the 
convergence of the optimal designs might be achieved by including analytically supplied 
gradient information to the optimizer, such as partial derivatives of the objective and 


constraint functions. 


Zt 











Table 5. Optimization Results for 6 Different Initial Design Values 


Starting Values (ho) for all 12 Design Variables 
230.9516 | 230.9296 | 233.8899 


Active Constraints 
Function Calis 











2. Solution Using Orthogonal Array Optimization Technique 


The same problem is solved again using the technique developed in this research. 
The optimization problem is changed slightly for implementation of this technique. Instead 
of keeping all of the constraints negative, the normalized stress, displacement, and lowest 
natural frequency are maintained between O and 1, where values of 1 denote active 
constraints. 

Therefore, the problem statement for the constrained optimization problem is given 


by the following objective function, F(h), and constraints, g(h): 














12 
minimize: F(h)= > bA,L (21) 
=] 
. |o 
subject to: g(i) =|—™™ <1 fori = 1,...,12 (22) 
yield 
aye (23) 
OW, 
Uu 
e(14)=—-- <1 (24) 
0.1<h, <4.0 fori = 1,...,12 (25) 
where: Uae = 0-12) 
Q.=7 Hz 


Syield = 78000 psi 
b; = 0.5 hj 
h; are the 12 design variables 


The constraints are calculated from FEM calculations 


29 





The first step is to select an appropriate orthogonal array to deal with the 12 design 
variables in this problem. Therefore, the Lj6(2!>) Orthogonal Array (Taguchi, 1987) as 


seen in Table 6 is used. 


Table 6. L46(2!5) Orthogonal Array 


L16(2+°) 
Parameter No. | 1 2 4 5 6 7 8 9 10 11 12 13 14 15 


No. of Runs 
1 l 1 1 ] 1 1 1 1 1 | 1 1 
2 1 1] 1 1] 1 J 1 2 2. 2 Zz 2@ 2 2 
3 I 1 |] 2. 2 Z2- 2 1 1 1 |] 2 2 2 2 
4 | 1 1 e+ Z 2 2 2 2 2 1 1 1 1 
5 l 2. 2 1 1 Zz 2 uf 2 2 | Ae 2 2 
6 1 2 2 1 1 2 2 a 2 1 1 2. 9 1 1 
7 1 2° oy 1 1 1 | 2. “2 Zz. 2 1 1] 
8 1 Zz 2 2D. 2 1 1] 2 2 1 1 1 | 2 2 
9 2 1 2 | ee! | 2 1 2 t 2 Lae 1 2 
10 ps 1 2 1 2 i 2 2 | 2 4 2 1 2 1 
11 2 1 2 2 4 ee. 1 2 1 2 2 | 2 1 
12 2 1 2 2. \4 2Z. I 2. ‘4 2 /I 1 2 t 2 
13 2 2 1 1 2 2 1 1 2 oe od 1 2 2 1 
14 2 2: ‘4 1 ‘2 Zz “| 2 1 1 2 2 I 2 
15 2 2 4 2 1 { 2 1 2 2 | 2 ft 1 2 
16 2 2. .'l 2 1 1 2 2 |i : 2 1 2 2 1 


The values of 1 in the orthogonal array will correspond to h = 1.00", and the values 
of 2 will correspond to h = 4.00". The FEM programs are then run for the 16 runs in the 
orthogonal array with the appropriate values of h for each design variable. Since there are 
only 12 design variables, the last 3 columns could be used for interactions, but will not be 
used in this example. The results from the running of the orthogonal array are shown in 


Table 7. 
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Table 7. Results from the L16(2!5) Orthogonal Array 


Norm. Norm. Norm. Norm. Norm. Norm. Norm. Norm. Norm. Norm. Norm. Norm. Norm. 
Stress. Stress. Stress. Stress. Stress. Stress. Stress. Stress. Stress. Stress. Stress. Stress. Norm. Min. 

Element Element Element Element Element Element Element Element Element Element Element Element Displ. Natural 
8 10 11 12 Freq. 
2.6902 2.0322 0.9896 1.3444 -0.0128 -0.0128 -0.0128 -0.0128 -1.3352 -0.9805 -1.9863 -2.6441 38.4095 2.2742 
0.1501 0.1827 0.2315 0.1737 -0.0546 -0.0546 -0.4163 -0.0116 -0.0305 -0.0631 -0.1081 -0.1774 1.8670 1.0978 
0.8151 0.8439 1.2901 0.0578 0.0279 0.0004 -0.0082 -0.3943 0.0210 0.1493 -2.2303 -0.0963 6.8973 2.9096 
0.6602 0.6581 0.1815 0.0086 -0.0083 -0.0200 -0.0126 -0.0516 -0.0454 -0.0196 -0.0963 -2.0859 4.0895 3.0482 
1.8663 0.0808 0.0368 0.8278 -0.0133 -0.0006 -0.0354 -2.2470 -0.6759 -0.0151 -0.0275 -0.8989 8.5306 2.9907 
1.14958 0.0560 0.0298 0.6499 -0.5278 -0.0088 -0.0270 -0.0270 -0.0258 -0.5912 -1.2226 -0.0573 3.4959 2.1182 
0.1583 0.0056 0.0155 0.0324 -0.0123 -0.5064 -0.0591] -0.0591 -0.5064 -0.0191 -0.0960 -0.1737 1.8519 1.4112 
2.1311 0.0896 0.0871 0.1003 -0.0330 -2.0618 0.4048 0.0066 -0.0330 -0.8608 -0.8804 -0.8828 8.7463 3.7037 
0.0855 1.8414 0.0450 1.0612 -0.0039 -0.2402 -0.0224 -1.4102 -0.0092 -0.2260 -0.0251 -0.7464 6.0302 2.5746 
0.0571 1.2030 0.0323 0.7029 -0.0030 -0.1607 -0.0054 -0.0054 -0.5441 -0.0235 -1.1767 -0.0553 3.2592 1.7746 
0.0039 0.1415 0.0090 0.0164 -0.7389 -0.0124 -0.0810 -0.0810 -0.0383 -0.8123 -0.0827 -0.1541 2.5589 1.9492 
0.1001 2.1992 0.0099 O.0115 -0.2700 -0.0054 1.1727 0.0192 -1.1891 -0.0517 -0.5254 -0.5843 6.9061 2.9060 
0.1426 0.0715 0.0543 0.0732 -0.0021 -0.0359 0.0294 0.0294 -0.0708 -0.0353 -0.0146 -0.0429 0.9212 1.1302 
0.0780 0.0460 1.1533 1.1723 -0.0008 -0.0008 -0.3639 -0.0064 -1.1679 -1.1506 -0.0426 -0.0822 8.0068 2.8112 
0.0165 -0.0297 -0.0297 0.0150 0.8846 -0.0161 -0.0326 0.1851 0.0033 4.4952 2.0479 
0.0243 -0.2669 0.0238 0.0222 -0.9740 -0.1114 7.8706 1.9429 








Weight | 














Run No. 
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From the FEM results for each run in the orthogonal array, the averaged Jacobian is 
computed (See Equation 3) which contains the average overall sensitivities of each of the 12 
design variables to the 14 constraints. The assumption is made that the sensitivities of the 
normalized stresses in each member due to other members is small, and therefore neglected, 
compared to the sensitivities of the normalized stresses due to itself. Therefore the averaged 
Jacobian matrix is simplified as shown in Table 8. 

The absolute maximum value in each column is then used to determine the amount 


each design variable, hj, should be changed because the absolute maximum value is the 


most constraining. It is denoted by 2 . Equation 6 will be used to calculate the 
j max 


value Ah;, where Ah; = Ax; . 

Table 9 shows the 1-D search for the hp=[4 444444444 4 4] starting point. 
The highlighted values are the maximum constraint values (1 is the target). For the case 
when P, = 2.2, the Maximum Normalized Displacement exceeds 1, therefore the distance in 
the search direction is bracketed and simple interpolation is used to get a value close to 1. 

The results for several different starting points are shown in Table 10. The case 
when hg = [11111 1 111 1 1] is initially in the infeasible domain, with a resulting 
value in the feasible domain, hj=[44444444 4 4 4 4], from which the procedure can be 
repeated. 

Figures 17-20 show the displacements of the frame for the 4 initial designs. The 
magnification of the displacements for these plots is 50x. Figures 21-23 show the 
displacements of the final design determined by using the orthogonal array for 3 of the initial 
conditions. 

Figures 24-26 show the first 3 mode shapes of the final design found from the initial 
design of ho=[4 44444444444]. The magnification of displacements is 10x. 
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Table 8. Average Overall Sensitivities from the L16(215) Orthogonal Array 
(a a a GR ae eC ee 

















0.14 
mt tT | Tt | TE tT 
AP. -0.03 
mt tT TT Td 
AP, 0.13 
m | | | | | dC ae 
AP, 0.25 
mt | tt tt | — 
an 





| 1.41] -1.09] -1.30] -1.13 | -1.43] -1.30] -1.03] -1.06] -2.06] -2.09] -1.35] -2.04 
Ah, 


= the sensitivity of the stress in element 1 to height 1 





< 
cf 
a 


= the sensitivity of the maximum displacement to height i 






* 








re = the sensitivity of the minimum natural frequency to height i 


f 


for i= 1 through 12 


35 








onal Array with hp=[44444444444 4] 

























Maximum Minimum Maximum Weight 

Stress Normalized Normalized 

Constraint Natural Displacement 
Frequenc 









0.6804 0.2898 537.96 









0.0802 0.7295 0.3784 471.93 
0.0903 0.7672 0.4579 430.50 
-0.1348 0.9144 0.8748 318.63 
-0.1573 0.9795 1.1195 285.49 
0.1455 0.9457 0.9874 301.80 











Figures 27-29 show the first 3 mode shapes of the final design found from the initial 
design of hop=[3 3 3 33 33 3 33 3 3]. The magnification of displacements is 10x. 

Figures 30-32 show the first 3 mode shapes of the final design found from the initial 
design of ho=[2 22222222222]. The magnification of displacements is 10x. 

As shown, the averaged Jacobian is a matrix of the average overall sensitivities 
which yields a search direction that, when coupled with a 1-D search technique, can 
determine a design which lowers the objective function and always remains in the feasible 


domain. The starting point value X%! is updated according to 
X47 = X1" + aS! (26) 
where X4 is the new vector of design variables, q is the iteration number, S is a search 


direction vector (here calculated by the averaged Jacobian matrix), and the scalar quantity o 


is the distance in the search direction, S. 
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Figure 17. Displacement for Initial Design, h=4.0, Max. Disp. = 0.0195" (Mag = 50x) 
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Figure 18. Displacement for Initial Design, h=3.0, Max. Disp. = 0.0595" (Mag = 50x) 
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Mesh of frame for 12 elements 
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Figure 19. Displacement for Initial Design, h=2.0, Max. Disp. 
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Figure 20. Displacement for Initial Design, h=1.0, Max. Disp. = 4.6091" (Mag = 50x) 
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Figure 21. Displacement for Final Design, ho=4.0, Max. Disp. = 0.1185" (Mag = 50x) 
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Figure 22. Displacement for Final Design, ho=3.0, Max. Disp. = 0.1206" (Mag = 50x) 
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Figure 23. Displacement for Final Design, hg=2.0, Max. Disp. = 0.1185" (Mag = 50x) 
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First Mode Shape, 7.40 Hz 
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Figure 24. Ist Mode Shape with ho=4.0 (Mag = 10x) 


Second Mode Shape, 27.4 Hz 
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Figure 25. 2nd Mode Shape with hp=4.0 (Mag = 10x) 


Third Mode Shape, 59.5 Hz 
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Figure 26. 3rd Mode Shape with ho=4.0 (Mag = 10x) 
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First Mode Shape, 7.56 Hz 
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Figure 27. Ist Mode Shape with hp=3.0 (Mag = 10x) 


Second Mode Shape, 25.44 Hz 
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Figure 28. 2nd Mode Shape with hp=3.0 (Mag = 10x) 


Third Mode Shape, 59.24 Hz 
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Figure 29. 3rd Mode Shape with ho=3.0 (Mag = 10x) 
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First Mode Shape, 7.70 Hz 
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Figure 30. 1st Mode Shape with ho=2 (Mag = 10) 
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Figure 31. 2nd Mode Shape with hg=2 (Mag = 10) 


Third Mode Shape, 60.08 Hz 





Figure 32. 3rd Mode Shape with ho=2, (Mag = 10) 
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Starting Values (ho) for all 12 Design Variables 


Active Constraints 
Function Calls 





In this example, a simple 1-D search algorithm is used which extrapolates out in the 
direction chosen until a constraint is violated. Once the constraint is violated, the solution is 


bounded and an interpolation procedure similar to the golden section method may be used. 
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3. Comparison of the two Optimization Techniques 


Table 11 shows a comparison of the convergence rates between the SQP 
optimization method and the Orthogonal Array Technique. It is important to notice that the 
SQP routine does not always approach the optimum from the feasible domain. However, 
the Orthogonal Array Technique always approaches the optimum solution from the feasible 


domain. This is an extremely important feature of the technique. 


ence Comparison 
Orthogonal Array 


Table 11. Conver 
SQP 




























Function Objective Max 
Evaluations Constraint 
Value <1 


23 301.801 9874 
22 285.57 1.005 
23 284.45 0.9878 


17 720 0.5872 
40 (17+23) 301.801 9874 


Function Objective Max 
Evaluations Constraint 
Value <0 


28 444.596 -0.5136 
43 326.061 -0.2585 
58 274.131 -0.554 

28 300.724 -0.2604 
44 283.249 -0.1643 
13 180 1.344 * 
351 230.875 0.0013 























ho = {1.0} 





13 35.871 * 
496 230.930 .0008 





* indicates infeasible design 


For the ho = {2.0} initial condition, the SQP converges to an optimum of 230.9 in 
351 function evaluations (with infeasible solutions at intermediate steps), and the 
Orthogonal Array method converges to 284.5 in only 23 function evaluations. The 
Orthogonal Array method converges to within 23% of the optimum design in only 6.5% of 
the number of function evaluations. This is a tremendous computational savings 


considering the FEM analyses required at each function evaluation. Further research is 
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being conducted to improve the convergence closer to the optimum, but studies so far have 


been inconclusive. 








V. CONCLUSION 


Improving design in engineering is a critical process in order for businesses to stay 
competitive and profitable. Therefore, optimization techniques are required to improve the 
design process and help the engineer to develop better products in the future. This thesis 
has described a new design optimization technique that couples the orthogonal array 
experimental procedure with finite element analysis in an effort to develop an economical 
and efficient way of optimizing engineering designs. This paper has also shown some 
examples of the procedure which produced very good results. However, the proposed 
design optimization technique needs to be further developed and refined for applications to 
various design problems. Results also show that the orthogonal array technique converged 
towards the optimum at a faster rate than the SQP method for the frame example. 
Additionally, the goal of reducing the number of FEM calculations was achieved as well as 
providing the engineer with insight into the effects of changing the design variables on the 
constraints. This sensitivity information may prove very useful when developing future 


designs of a similar nature. 
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VI. RECOMMENDATIONS 


Further coupling of the orthogonal array technique with one-dimensional search 
techniques will provide for a more robust optimization algorithm and should be the topic of 
future research. In the 1-D search technique a possible technique might be to vary each 
value of the target value individually, in an effort to further reduce the objective function 
with a single application of the orthogonal array. 

Also, additional studies should look at the problem of refining the solution further 
towards the actual optimum. One possible method of doing this might be to incorporate a 
second orthogonal array around a refined domain. 

Studies should be conducted on additional highly non-linear problems to further 
evaluate the performance of this technique. 

The cases shown in this thesis, show designs that mostly start from the feasible 
domain. Additional research should be conducted for cases of an initially infeasible design. 

Additionally, utilizing all of the average overall sensitivities and not just the diagonal 
terms in the first n rows of the.averaged Jacobian should be examined to enhance the 
procedure developed in this research. 

Finally, research conducted by Yurkovich (1994) suggests a similar approach to the 
optimization of structures with the utilization of Taguchi Methods. However, Yurkovich 
used Taguchi Methods to select the critical design variables and then implimented a 
conventional optimization technique. 

From Yurkovich's research, the incorporation of interactions into the procedure 
developed in this thesis might prove very useful in improving the optimization technique. 
Additionally, orthogonal arrays can be further utilized to reduce the number of design 


variables and interactions that must be analyzed with the finite element method. 


47 





48 








LIST OF REFERENCES 


Ealey, Lance A. 1988. Quality by Design. Dearborn, MI: ASI Press. 

Grace, Andrew. 1992. Optimization Toolbox User's Guide Natick, MA: The MathWorks, 
Inc. 

Kwon, Young W. and Hyochoong Bang. 1997. The Finite Element Method using 
MATLAB. Boca Raton, FL: CRC Publications. 

Lang, Serge. 1968. Calculus. New York: Addison-Wesley. 

Peace, Glen S. 1993. Taguchi Methods: A Hands-On Approach To Quality Engineering. 
Reading, MA: Addison-Wesley. 

Phadke, Madhav S. 1989. Quality Engineering using Robust Design. Englewood Cliffs, 
NJ: Prentice-Hall. 

Ryan, Nancy E., ed. 1988. Taguchi Methods and QFD. Dearborn, MI: ASI Press. 

Taguchi, G. and Konishi, S. 1987. Taguchi Methods: Orthogonal Arrays and Linear 
Graphs. Dearborn, MI: ASI Press. 

Vanderplaats, Garret N. 1984. Numerical Optimization Techniques for Engineering Design: 
With Applications. New York: McGraw-Hill. 

Yurkovich, Rudy. 1994. The use of Taguchi Techniques with the ASTROS Code for 
Optimum Wing Structural Design. AIAA Paper 94-1484-CP, 1994 Structures, 


Structural Dynamics, and Materials Conference. 


49 


10. 


INITIAL DISTRIBUTION LIST 


Defense Technical Information Center 
8725 John J. Kingman Rd., STE 0944 
Ft. Belvoir, VA 22060-6218 


Dudley Knox Library 
Naval Postgraduate School 
411 Dyer Rd. 

Monterey, CA 93943-5101 


Chairman, Code ME/Mc 

Department of Mechanical Engineering 
Naval Postgraduate School 

Monterey, CA 93943 


Professor Young W. Kwon, ME/Kw 
Department of Mechanical Engineering 
Naval Postgraduate School 

Monterey, CA 93943 


Naval Mechanical Engineering, Code 34 
Naval Postgraduate School 
Monterey, CA 93943 


COL (Ret) Stephen H. Young 
3620 Mountain View Dr. 
Boise, ID 83704 


Stuart H. Young 
8621 Rocky Gap Ct. 
Lorton, VA 22079 


Dr. Philip J. Emmerman 

US Army Research Laboratory 
2800 Powder Mill Rd. 
Adelphi, MD 20783 


Dr. Thomas K. Mills 

US Army Research Laboratory 
2800 Powder Mill Rd. 
Adelphi, MD 20783 


Barbara Broome 

US Army Research Laboratory 
2800 Powder Mill Rd. 
Adelphi, MD 20783 


51 


No. Copies 
) 





